function variational_properties()
% VARIATIONAL_PROPERTIES Hermitian矩阵特征值的变分性质
% 
% 展示Hermitian矩阵特征值的变分表示和相关性质

fprintf('=== Hermitian矩阵特征值的变分性质 ===\n\n');

%% 1. Rayleigh商基础
fprintf('1. Rayleigh商基础\n');
fprintf('目标：理解Rayleigh商的定义和性质\n\n');

% 对称矩阵
A1 = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('对称矩阵 A1 =\n');
disp(A1);

% 计算精确特征值
eigenvals1 = sort(eig(A1), 'descend');
fprintf('特征值（降序）: [%.6f, %.6f, %.6f]\n', eigenvals1);

% 测试不同向量的Rayleigh商
test_vectors = {
    [1; 0; 0],      % 标准基向量
    [0; 1; 0],
    [0; 0; 1],
    [1; 1; 1]/sqrt(3),  % 归一化向量
    [1; -1; 1]/sqrt(3),
    [2; 1; -1]/sqrt(6)
};

fprintf('\nRayleigh商测试:\n');
fprintf('向量                    Rayleigh商    与特征值的关系\n');
fprintf('----                    ----------    --------------\n');

for i = 1:length(test_vectors)
    x = test_vectors{i};
    x = x / norm(x);  % 归一化
    
    % 计算Rayleigh商
    rayleigh_quotient = (x' * A1 * x) / (x' * x);
    
    % 找到最接近的特征值
    [~, closest_idx] = min(abs(eigenvals1 - rayleigh_quotient));
    
    fprintf('[%4.2f,%4.2f,%4.2f]    %10.6f    接近λ_%d\n', ...
            x, rayleigh_quotient, closest_idx);
end

fprintf('\n观察：Rayleigh商总是在最小和最大特征值之间\n');
fprintf('范围: [%.6f, %.6f]\n', eigenvals1(end), eigenvals1(1));

%% 2. 变分原理：最大最小定理
fprintf('\n2. 变分原理：最大最小定理\n');
fprintf('目标：验证特征值的变分表示\n\n');

% 最大特征值的变分表示：λ₁ = max_{||x||=1} x'Ax
fprintf('最大特征值的变分表示：λ₁ = max_{||x||=1} x''Ax\n');

% 随机搜索验证
n_trials = 1000;
rayleigh_values = zeros(n_trials, 1);

for i = 1:n_trials
    x_random = randn(3, 1);
    x_random = x_random / norm(x_random);
    rayleigh_values(i) = x_random' * A1 * x_random;
end

max_rayleigh = max(rayleigh_values);
min_rayleigh = min(rayleigh_values);
mean_rayleigh = mean(rayleigh_values);

fprintf('随机搜索结果（%d次试验）:\n', n_trials);
fprintf('最大Rayleigh商: %.6f (理论值: %.6f)\n', max_rayleigh, eigenvals1(1));
fprintf('最小Rayleigh商: %.6f (理论值: %.6f)\n', min_rayleigh, eigenvals1(end));
fprintf('平均Rayleigh商: %.6f\n', mean_rayleigh);

% 最小特征值的变分表示：λₙ = min_{||x||=1} x'Ax
fprintf('\n最小特征值通过最小化Rayleigh商得到\n');
fprintf('误差：最大值 %.2e，最小值 %.2e\n', ...
        abs(max_rayleigh - eigenvals1(1)), abs(min_rayleigh - eigenvals1(end)));

%% 3. Courant-Fischer定理
fprintf('\n3. Courant-Fischer定理\n');
fprintf('目标：展示中间特征值的变分表示\n\n');

fprintf('Courant-Fischer定理：第k个特征值可以表示为\n');
fprintf('λₖ = min_{dim(S)=k} max_{x∈S,||x||=1} x''Ax\n');
fprintf('或者\n');
fprintf('λₖ = max_{dim(S)=n-k+1} min_{x∈S,||x||=1} x''Ax\n\n');

% 计算第2个特征值（中间特征值）
% 方法：在与第1个特征向量正交的子空间中最大化Rayleigh商
[V1, D1] = eig(A1);
[eigenvals1_sorted, sort_idx] = sort(diag(D1), 'descend');
V1_sorted = V1(:, sort_idx);

v1 = V1_sorted(:, 1);  % 第1个特征向量
fprintf('第1个特征向量 v1 = [%.4f, %.4f, %.4f]''\n', v1);

% 在v1的正交补空间中搜索
n_trials_2 = 1000;
rayleigh_orthogonal = zeros(n_trials_2, 1);

for i = 1:n_trials_2
    % 生成随机向量
    x_rand = randn(3, 1);
    
    % 投影到v1的正交补空间
    x_orth = x_rand - (x_rand' * v1) * v1;
    
    if norm(x_orth) > 1e-10  % 避免零向量
        x_orth = x_orth / norm(x_orth);
        rayleigh_orthogonal(i) = x_orth' * A1 * x_orth;
    end
end

% 移除零值
rayleigh_orthogonal = rayleigh_orthogonal(rayleigh_orthogonal ~= 0);
max_rayleigh_orth = max(rayleigh_orthogonal);

fprintf('在v1正交补空间中的最大Rayleigh商: %.6f\n', max_rayleigh_orth);
fprintf('第2个特征值（理论值）: %.6f\n', eigenvals1_sorted(2));
fprintf('误差: %.2e\n', abs(max_rayleigh_orth - eigenvals1_sorted(2)));

%% 4. 特征值的单调性
fprintf('\n4. 特征值的单调性\n');
fprintf('目标：展示矩阵扰动对特征值的影响\n\n');

% 原始矩阵
A4 = [3, 1; 1, 2];
eigenvals_orig = sort(eig(A4), 'descend');
fprintf('原始矩阵特征值: [%.6f, %.6f]\n', eigenvals_orig);

% 添加正定扰动
perturbations = [0.1, 0.5, 1.0, 2.0];
fprintf('\n扰动量    新特征值1    新特征值2    变化量1    变化量2\n');
fprintf('------    ---------    ---------    -------    -------\n');

for pert = perturbations
    % 正定扰动矩阵
    E = pert * eye(2);
    A_perturbed = A4 + E;
    
    eigenvals_pert = sort(eig(A_perturbed), 'descend');
    changes = eigenvals_pert - eigenvals_orig;
    
    fprintf('%6.1f    %9.6f    %9.6f    %7.4f    %7.4f\n', ...
            pert, eigenvals_pert, changes);
end

fprintf('\n观察：正定扰动使所有特征值增加（单调性）\n');

% 负定扰动
fprintf('\n负定扰动的影响:\n');
fprintf('扰动量    新特征值1    新特征值2    变化量1    变化量2\n');
fprintf('------    ---------    ---------    -------    -------\n');

for pert = perturbations
    E = -pert * eye(2);
    A_perturbed = A4 + E;
    
    eigenvals_pert = sort(eig(A_perturbed), 'descend');
    changes = eigenvals_pert - eigenvals_orig;
    
    fprintf('%6.1f    %9.6f    %9.6f    %7.4f    %7.4f\n', ...
            -pert, eigenvals_pert, changes);
end

%% 5. 交错定理（Interlacing Theorem）
fprintf('\n5. 交错定理\n');
fprintf('目标：展示子矩阵特征值与原矩阵特征值的关系\n\n');

% 4×4对称矩阵
A5 = [4, 1, 0, 1; 1, 3, 1, 0; 0, 1, 2, 1; 1, 0, 1, 3];
fprintf('4×4对称矩阵 A5 =\n');
disp(A5);

eigenvals_full = sort(eig(A5), 'descend');
fprintf('完整矩阵特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals_full);

% 提取3×3主子矩阵
A5_sub = A5(1:3, 1:3);
eigenvals_sub = sort(eig(A5_sub), 'descend');
fprintf('3×3主子矩阵特征值: [%.4f, %.4f, %.4f]\n', eigenvals_sub);

% 验证交错性质：λ₁ ≥ μ₁ ≥ λ₂ ≥ μ₂ ≥ λ₃ ≥ μ₃ ≥ λ₄
fprintf('\n交错性质验证:\n');
fprintf('λ₁ = %.4f ≥ μ₁ = %.4f: %s\n', eigenvals_full(1), eigenvals_sub(1), ...
        mat2str(eigenvals_full(1) >= eigenvals_sub(1)));
fprintf('μ₁ = %.4f ≥ λ₂ = %.4f: %s\n', eigenvals_sub(1), eigenvals_full(2), ...
        mat2str(eigenvals_sub(1) >= eigenvals_full(2)));
fprintf('λ₂ = %.4f ≥ μ₂ = %.4f: %s\n', eigenvals_full(2), eigenvals_sub(2), ...
        mat2str(eigenvals_full(2) >= eigenvals_sub(2)));
fprintf('μ₂ = %.4f ≥ λ₃ = %.4f: %s\n', eigenvals_sub(2), eigenvals_full(3), ...
        mat2str(eigenvals_sub(2) >= eigenvals_full(3)));
fprintf('λ₃ = %.4f ≥ μ₃ = %.4f: %s\n', eigenvals_full(3), eigenvals_sub(3), ...
        mat2str(eigenvals_full(3) >= eigenvals_sub(3)));
fprintf('μ₃ = %.4f ≥ λ₄ = %.4f: %s\n', eigenvals_sub(3), eigenvals_full(4), ...
        mat2str(eigenvals_sub(3) >= eigenvals_full(4)));

%% 6. 惯性定理（Sylvester's Law of Inertia）
fprintf('\n6. 惯性定理\n');
fprintf('目标：展示矩阵的惯性在合同变换下不变\n\n');

% 不定矩阵（有正负特征值）
A6 = [2, 1, 0; 1, 0, 1; 0, 1, -1];
fprintf('不定矩阵 A6 =\n');
disp(A6);

eigenvals6 = eig(A6);
n_positive = sum(eigenvals6 > 1e-10);
n_negative = sum(eigenvals6 < -1e-10);
n_zero = sum(abs(eigenvals6) <= 1e-10);

fprintf('特征值: [%.4f, %.4f, %.4f]\n', sort(eigenvals6, 'descend'));
fprintf('惯性：正特征值 %d 个，负特征值 %d 个，零特征值 %d 个\n', ...
        n_positive, n_negative, n_zero);

% 合同变换
P6 = [1, 2, 1; 0, 1, 1; 1, 0, 2];
A6_congruent = P6' * A6 * P6;
fprintf('\n合同变换后的矩阵 P''*A6*P =\n');
disp(A6_congruent);

eigenvals6_cong = eig(A6_congruent);
n_positive_cong = sum(eigenvals6_cong > 1e-10);
n_negative_cong = sum(eigenvals6_cong < -1e-10);
n_zero_cong = sum(abs(eigenvals6_cong) <= 1e-10);

fprintf('合同变换后特征值: [%.4f, %.4f, %.4f]\n', sort(eigenvals6_cong, 'descend'));
fprintf('合同变换后惯性：正特征值 %d 个，负特征值 %d 个，零特征值 %d 个\n', ...
        n_positive_cong, n_negative_cong, n_zero_cong);

fprintf('惯性保持不变: %s\n', mat2str(n_positive == n_positive_cong && ...
                                    n_negative == n_negative_cong && ...
                                    n_zero == n_zero_cong));

%% 7. 特征值的包含域
fprintf('\n7. 特征值的包含域\n');
fprintf('目标：使用Gershgorin圆盘定理估计特征值范围\n\n');

A7 = [5, 1, 0.5; 0.2, 3, 0.8; 0.1, 0.3, 2];
fprintf('矩阵 A7 =\n');
disp(A7);

eigenvals7 = sort(eig(A7), 'descend');
fprintf('精确特征值: [%.4f, %.4f, %.4f]\n', eigenvals7);

% Gershgorin圆盘
fprintf('\nGershgorin圆盘:\n');
for i = 1:size(A7, 1)
    center = A7(i, i);
    radius = sum(abs(A7(i, :))) - abs(A7(i, i));
    fprintf('圆盘 %d: 中心 = %.4f, 半径 = %.4f, 区间 = [%.4f, %.4f]\n', ...
            i, center, radius, center - radius, center + radius);
end

% 验证特征值是否在圆盘内
fprintf('\n特征值包含验证:\n');
for i = 1:length(eigenvals7)
    lambda = eigenvals7(i);
    contained = false;
    for j = 1:size(A7, 1)
        center = A7(j, j);
        radius = sum(abs(A7(j, :))) - abs(A7(j, j));
        if abs(lambda - center) <= radius + 1e-10
            contained = true;
            fprintf('λ_%d = %.4f 包含在圆盘 %d 中\n', i, lambda, j);
            break;
        end
    end
    if ~contained
        fprintf('λ_%d = %.4f 不在任何圆盘中（不应该发生）\n', i, lambda);
    end
end

%% 8. 变分原理的应用：特征值优化
fprintf('\n8. 变分原理的应用：特征值优化\n');
fprintf('目标：通过变分原理求解约束优化问题\n\n');

% 问题：在约束 x'Bx = 1 下最大化 x'Ax
A8 = [3, 1; 1, 2];
B8 = [2, 0.5; 0.5, 1];

fprintf('目标矩阵 A8 =\n'); disp(A8);
fprintf('约束矩阵 B8 =\n'); disp(B8);

% 这等价于广义特征值问题 A8*x = λ*B8*x
[V8, D8] = eig(A8, B8);
eigenvals8 = diag(D8);
[max_eigenval, max_idx] = max(eigenvals8);
optimal_x = V8(:, max_idx);

fprintf('广义特征值: [%.6f, %.6f]\n', eigenvals8);
fprintf('最大值: %.6f\n', max_eigenval);
fprintf('最优解: x = [%.4f, %.4f]''\n', optimal_x);

% 验证约束和目标函数值
constraint_value = optimal_x' * B8 * optimal_x;
objective_value = optimal_x' * A8 * optimal_x;

fprintf('约束验证 x''Bx = %.6f (应该接近1)\n', constraint_value);
fprintf('目标函数值 x''Ax = %.6f\n', objective_value);
fprintf('比值 (x''Ax)/(x''Bx) = %.6f (应该等于最大广义特征值)\n', ...
        objective_value / constraint_value);

%% 总结
fprintf('\n=== 变分性质总结 ===\n');
fprintf('1. Rayleigh商提供特征值的变分表示\n');
fprintf('2. 最大最小定理给出特征值的优化特征\n');
fprintf('3. Courant-Fischer定理处理中间特征值\n');
fprintf('4. 特征值关于矩阵扰动具有单调性\n');
fprintf('5. 交错定理描述子矩阵与原矩阵特征值的关系\n');
fprintf('6. 惯性定理说明惯性在合同变换下不变\n');
fprintf('7. Gershgorin定理提供特征值的包含域\n');
fprintf('8. 变分原理可用于约束优化问题\n\n');

fprintf('这些性质在特征值计算和分析中具有重要意义：\n');
fprintf('• 提供特征值的理论界限\n');
fprintf('• 指导数值算法的设计\n');
fprintf('• 分析扰动对特征值的影响\n');
fprintf('• 解决约束优化问题\n\n');

end